 module simulate
      use precisions
      use types_generic
      use routines_specific
    contains

    subroutine dosim(grids, params, tech, sim, valpol, profiles, run)
        implicit none
        
        !! Argument
        type(typeGrids), intent(in) :: grids
        type(typeParams), intent(in)  :: params      
        type(typeTech), intent(in) :: tech
        type(typeSim), intent(in) :: sim   
        type(typeValPol), intent(in) :: valpol
        type(typeProfiles), intent(out) :: profiles
        integer, intent(in) :: run
        
        !! For program
        type(typeInfoForSim) :: infoForSim
        type(typeInfoForMax) :: infoForMax
        real (kind=rk) :: valstar, optPol(6)
        real (kind=rk) :: annuityInc, annnuityPot , weights(grids%maxT, sim%numSims) 
        type(typeForOptRout) :: ForOptRoute
        
        integer ::  scratch 
        character (130) :: path_in,  path_out, outfile
 
        !! Indexes
        integer :: ixT, ixS

   
        do ixS = 1, sim%numSims
 

             annuityInc = 0.0_rk
             infoForMax%ixA = -999 ! This is used only in printing in getMax. Set in solvVSL. Don't set here
             infoForMax%solveT_SimF = .false. !Set here instead of outside loop because of parallisationz
                            
            do ixT = 1, grids%maxT              
                infoForMax%today = ixT
                infoForMax%tomorrow = ixT + 1

               
                !--------------------!
                ! Set asset variables                
                !--------------------!
                if (ixT.eq.1) then
                    profiles%assStart(ixT, ixS) = sim%startA                  
                else   
                    profiles%assStart(ixT, ixS) = profiles%assEnd(ixT-1, ixS)  * (1+tech%r)
                end if
                

                
                infoForSim%startA= profiles%assStart(ixT, ixS)
                

                !-----------------------------!
                ! Set annuity income variable
                !-----------------------------!
                if (ixT .lt. (tech%retireAge )) then
                        infoForMax%annY = 0.0_rk
                elseif (ixT .ge. (tech%retireAge)) then ! In the last year of work (the year in which one can annuitize) income goes to whatever the new annuity stream will be
                          infoForMax%annY = annuityInc 
                    
                end if        
                infoForMax%Y =  infoForMax%annY                            

                if (ixT .ge. (tech%retireAge)) then ! In the last year of work (the year in which one can annuitize) income goes to whatever the new annuity stream will be
                    infoForMax%cashOnHand = profiles%assStart(ixT, ixS) + infoForMax%Y  + tech%SSlevel
                else
                    infoForMax%cashOnHand = profiles%assStart(ixT, ixS) + infoForMax%Y                  
                end if
                
                infoForMax%Y1 = infoForMax%annY 

                call findingridr(grids%A(:,ixT),  grids%maxA, infoForMax%cashOnHand, scratch, profiles%assLocStart(ixT, ixS))
                
                
                call getMax(infoForMax, grids, valpol, params, tech, optPol, valstar)

                
                profiles%cons(ixT, ixS) = optPol(1)
                profiles%val(ixT, ixS) = valstar
                if (ixT .eq. (tech%retireAge - 1)) then
                    annnuityPot = optPol(2) * (infoForMax%cashOnHand)
                    annuityInc =  annnuityPot * tech%annRates(infoForMax%today) 
                    profiles%annProp(:, ixS) = annnuityPot / (infoForMax%cashOnHand - profiles%cons(ixT, ixS))
                   
                    if (printingGlob.eq.2) write(*,*) ixS, ixT, profiles%annProp(ixT, ixS)
                   
                    profiles%assEnd(ixT, ixS) = (infoForMax%cashOnHand - annnuityPot ) - profiles%cons(ixT, ixS) 
                     if (profiles%cons(ixT, ixS) .gt. 0.000_rk) then
                        write(*,*) 'cons should be zero and isnt in simulate'
                        stop
                    end if
                
                        
                  
                    if ((infoForMax%cashOnHand   - profiles%cons(ixT, ixS))   .lt.0.001) then
                        weights(:, ixS) = 0.0_rk
                       if (profiles%cons(ixT, ixS) .gt. 0.000_rk) then
                        write(*,*) 'weight should not be zero and isnt in simulate'
                        stop
                    else 
                        weights(:, ixS) = 1.0_rk
                    end if
                    end if
                
                else
                    profiles%assEnd(ixT, ixS) = (infoForMax%cashOnHand )  - profiles%cons(ixT, ixS)           
                end if
                
                
               if (ixT.lt. grids%maxT) then
                    call findingridr(grids%A(:,ixT+1),  grids%maxA, profiles%assEnd(ixT, ixS), scratch, profiles%assLocEnd(ixT, ixS))
               else
                    profiles%assLocEnd(ixT, ixS) = 0
               end if
               

                
                profiles%Y(ixT, ixS) = infoForMax%Y
                profiles%annY(ixT, ixS) = annuityInc
                if (tech%annOn) then
                    profiles%annW(ixT, ixS) = annuityInc / tech%annRates(infoForMax%today)
                else
                    profiles%annW(ixT, ixS) = 0.0_rk
                end if                    
                profiles%totW(ixT, ixS) = profiles%assEnd(ixT, ixS) +  profiles%annW(ixT, ixS) 
                profiles%totLiqW(ixT, ixS) = profiles%assEnd(ixT, ixS)

            
            profiles%assLocStartReal(ixT, ixS) = real(profiles%assLocStart(ixT, ixS), rk)
            profiles%assLocEndReal(ixT, ixS) = real(profiles%assLocEnd(ixT, ixS), rk)

            end do
        end do
  
        






!write(*,*) 
write(*,'(6f12.3)') params%beta, params%gamma, profiles%annProp(45-19, 1), tech%annRates(1), profiles%Y(51, 1), profiles%cons(2, 1) !, advance = 'no'




        call getpaths(path_in, path_out, run)                    
        outfile = trim(path_out) // 'results.out'               
        open (unit=250, file=outfile, recl=25000, action='write', access = 'append')
        write (250, '(6f12.4)',  advance = 'no') params%beta, params%gamma, profiles%annProp(45-19, 1), tech%annRates(1), tech%annRateSubjForPrinting(1), profiles%Y(51, 1)
        close(250)
!write(*,*) 



    end subroutine dosim
    
    
    end module simulate
